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Abstract 

We investigate the nonequilibrium behavior of dense, plastic-crystalline suspensions of mildly 
anisotropic colloidal hard dumbbells under the action of an oscillatory shear field by employing 
Brownian dynamics computer simulations. In particular, we extend previous investigations, where 
we uncovered novel nonequilibrium phase transitions, to other aspect ratios and to a larger nonequi¬ 
librium parameter space, that is, a wider range of strains and shear frequencies. We compare and 
discuss selected results in the context of novel scattering and rheological experiments. Both simu¬ 
lations and experiments demonstrate that the previously found transitions from the plastic crystal 
phase with increasing shear strain also occur at other aspect ratios. We explore the transition 
behavior in the strain-frequency phase and summarize it in a nonequilibrium phase diagram. Ad¬ 
ditionally, the experimental rheology results hint at a slowing down of the colloidal dynamics with 
higher aspect ratio. 
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I. INTRODUCTION 


The equilibrium and nonequilibrium behavior of colloidal hard-sphere suspensions has 
been investigated extensively CHI! as it constitutes the simplest model system to gain a 
better understanding of the structure and phase behavior of colloidal suspensions in some 
of their fundamental aspects. Notwithstanding that the particles interact isotropically by 
excluded volume only, hard sphere colloids show a nontrivial phase diagram already in 
equilibrium and glassy behavior above a certain volume fraction muu- Regarding their 
nonequilibrium behavior, hard sphere suspensions subjected to shear have been of particular 
interest, since shear is one of the most common external fields Here, many interest¬ 

ing forms of order-disorder transitions have been observed in experiments and simulations 
under various shear conditions um- These include steady shear (constant rate) [H], oscil¬ 
latory shear (fixed strain amplitude and frequency) (Tj |8| , and oscillatory fixed shear rate 
protocols [9|. 


In particular, various nonequilibrium states have been identified in spherical colloids 
under shear and have been well characterized up to now using a combination of scattering 
methods [8], optical techniques mm, simulations [7] and theoretical models m- Similar 
order-disorder transitions of colloids under shear have been reported for clusters of soft 
spheres 


and highly charged spheres Most recently, [Besseling et al.\ combined 

confocal microscopy and Brownian dynamics (BD) simulations to explore the nonequilibrium 
behavior of colloidal hard spheres under oscillatory shear [7]. Here, for small strains a face- 
centered cubic (FCC) twin is found to be stable corroborating with the classical results [8|. 
At high strains the predominant structure is found to be registered sliding of hexagonal close 
packed (HOP) layers. Moreover, a dense direction of the HOP layer prefers to be parallel to 
the velocity at high strain amplitudes, whereas at low strains a dense direction is parallel 
to the flow direction. This behavior results in a 30° turn in the scattering pattern ng. 


The corresponding diffraction patterns show three-fold symmetries. [Besseling et al. 


have calculated an extensive nonequilibrium state diagram for hard spheres under oscillatory 
strain, which categorizes further high-strain structures. 


In general, however, colloids are anisotropic and as such there is a rising interest in the 
behavior of suspensions of non-spherical particles for fundamental understanding pTl - fTS] or 
applications, such as constituents for novel materials ng or photonics |2nH22|- Desirable 
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on a fundamental level are hard particles with slight anisotropy that weakly perturb the 
isotropic interactions of spherical reference systems. One popular experimental realization 
is a system of steeply repulsive ’dumbbells’, that is a colloidal dimer made up by two fused 
equally-sized hard spheres [23l [2lj . The equilibrium phase diagram and the stability and 
nucleation processes of (plastic) crystal phases of hard dumbbells have been mapped out 
comprehensively by means of Monte Carlo (MC) simulations |25H28j . Plastic crystal phases 
are characterized by a crystalline center of mass order and a lack of long-range order of 
the particles’ orientations [2SHHI], in which the latter is the distinguishing attribute in 
comparison to a fully ordered crystal. For the convenience of the reader the relevant parts 
of the phase diagram are replotted in Fig. 

Out of equilibrium, the impact of the weak anisotropy and its accompanying translational- 
rotational coupling ra on the rheology and the structure of dispersions under shear has been 
of keen interest. This aspect has been investigated by means of mode-coupling theory m 
[HH], neutron scattering [2H] and rheological experiments [311 EH]- A few years ago, we 
have presented a very neat experimental realization [HH] of monodisperse hard dumbbells 
(aspect ratios 0.24 and 0.30), which matches the plastic crystal phase boundaries well in the 
estimated phase diagram (Fig. [^. Employing this well-defined experimental model system, 
we have, in fact, very recently shown by a combination of rheology-scattering (rheo-SANS) 
and Brownian dynamics (BD) simulations that hard plastic-crystalline dumbbells undergo 
nonequilibrium transitions under oscillatory shear similar than hard spheres but the nature 
of the transition differs strikingly IS]. In particular, we observed a finite orientational 
correlation arising on increasing strain amplitudes [37] and a more vigorous transition with 
large implications on rheology and the yielding behavior. 

In this work, we follow up on our previous work and extend our computational work to 
other aspect ratios and a wider frequency range. Our focus is thus the region of small to 
moderate elongations (aspect ratios L* < 0.4) and high volume fractions where the plastic 
crystal (PC) phase predominates [3H|- Furthermore, kinetics properties of the suspensions 
are investigated by our simulations. We also present new rheology experiments on a second 
aspect ratio (0.3). We note that the aspect ratio of 0.3 represents an anisotropy very close 
to that of a nitrogen molecule, that also features a plastic crystal phase [39l 00], referred to 
as /5-phase, so one could consider our system as colloidal nitrogen. 
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FIG. 1. (Color online) Hard dumbbell phase diagram in the volume fraction (j) to aspect ratio L* 
plane [27]. The state points for dumbbells A (0 = 0.55, L* = 0.08), B (0 = 0.55, L* = 0.24), and 
C (0 = 0.55, L* = 0.3), considered in this work are marked by symbols (■). The state point S (•) 
denotes an almost hard-sphere-like reference system (0 = 0.55, L* = 0.02). For better orientation, 
the freezing (^Seit) close packed (0Qp) packing densities of the HS system 

are indicated at the vertical axis. 

II. METHODS 

A. Brownian Dynamics 

The BD simulations are carried out using Ermak ’s HU method for interacting particles 
in solution with an additional term to account for the oscillatory shear force. We neglect 
hydrodynamic interactions among the particles as the systems under consideration are dis¬ 
turbed in a regime where the driving forces are in the order of the viscous forces, i.e., the 
Peclet numbers are small, see also related work [Tj- The dumbbell particles are represented 
by a two-segment Yukawa model, in which the particles interact via two spherical beads con¬ 
straint at the constant distance L with a steep Yukawa potential (decay length = 0.05cr) 
defined by 

V(r) = €—exp{—K(r — a)} , (1) 

r 

where e = ksT and a set the energy and length scales, respectively, and r denotes the center- 
to-center distance between two beads. The parameter k tunes the softness of the interaction 
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and is chosen to maintain computational performance and resemble hard particle behavior. 
We employ a forward Euler scheme in order to integrate the equations of motion [HI |12] . 
The time scale is set to the Brownian time r = <j^/Dq of a single bead of diameter a and 
diffusivity Dq = kBT{3TTr]sa)~^ with a solvent viscosity rjs. The parallel and perpendicular 
center of mass (COM) coordinates are updated according to 
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The shear flow only affects the COM transport in a:-direction through the last term in 
equation (Eq. 0) The directors are updated following 
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where T” is the total torque exerted on particle i at time t = nAt. The single-particle 
diffusion coefficients parallel D\\ , perpendicular D± to the long axis and Dr for the rotation 
about the short axis depend on the particle geometry, their respective values are listed below 
(Table |I]). The torque is comprised of the inter-particle and background-flow contributions, 
via 


T,(t) = Tr(t) - ^ {u,(t) X r(() ■ u,(i)}. 


( 6 ) 


We impose a time-dependent linear shear flow in the a:-direction, such that the shear 
gradient is parallel to By and the vorticity is in 2 ;-direction. As the flow velocity v(?/) = 
Vx{y)^x vanishes at y = 0 and depends linearly on y, the velocity-vorticity plane is the 
plane of (spatially) constant flow velocity. Hence, the sheared system retains symmetry in 
this particular plane, cf. Fig. Mathematically, the shear flow is thus described by the 
rate-of-strain tensor 
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FIG. 2. (Color online) Sketch of the imposed shear flow and the used coordinate system where 
the velocity is in x, the velocity gradient in and the vorticity is along the 2 ;-direction. The 
dumbbells’ beads diameter is a and the aspect ratio (elongation) is defined as L* = L/a, where L 
is the center-to-center distance of both beads. In the simulations, we apply Lees-Edwards boundary 
conditions 031. 


which may be written as a sum of symmetric and anti-symmetric tensors, i.e., the shear and 
vorticity tensors 


E{t) = T{t) + n{t) 
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(9) 


The symmetric part T describes the pure elongational flow, while the anti-symmetric part 
e describes the vorticity of the flow. The vorticity is spatially uniform and quantifies the 
solvent contribution to the rotation of the particles. The time-dependent strain with the 
dimensionless amplitude imposed by the flow field on the suspension is thus given by 


lit) = 7maxSin(a;t). 

Hence, the linear solvent velocity profile is 

Vxiy,t) = iit)y. 


( 10 ) 


( 11 ) 


In order to compare the driving force to the intrinsic viscous forces, we define the Peclet 
numbers 

P® — /r^o ~ and 

f Imax./Dr 1 
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L* 

D\\lDl 

D±/Dl 

Dr/D^ 

0.02 (S) 

0.99 

0.99 

0.97 

0.10 (C) 

0.97 

0.95 

0.85 

0.24 (A) 

0.93 

0.89 

0.69 

0.30 (B) 

0.91 

0.87 

0.63 


TABLE 1. Single particle diffusive properties used in the simulations, obtained by SHM calculations. 


where the latter is a definition respecting the time scale set by the rotational Brownian 
motion. The maximum shear rate y^ax = 7max/'?' sets the time scale of the driving force 
exerted by the shear flow. 

In order to compare the results of our simulations, the time scales set by the diffusion 
constants and sizes of the particles are important. For dumbbells with L* 0.24 the 
parallel (Dy = 0.93 T)q), perpendicular = 0.89T)q) and rotational = 0.69D^ = 
0.69{3Dq)/{2Rh)‘^) diffusivities have been obtained by the shell bead model method (SHM) 
and matched to the experimental data described in our previous work [SH! • In Table |I] the 
reader may find the exact parameters used for the respective systems. Thus, the rotational 
Brownian time scale is Tr = 9.Q9{2 RhY/{3Dq) in this case. In length and time units of 2Rh 
and r the translational Peclet number is defined as Pe = ^/7max(2i?ir)^/T>f, where Rh is 
the hydrodynamic radius of a sphere in the experimental frame of reference. 

We simulate (A^ = 864) dumbbell particles subjected to Lees-Edwards [IH] periodic 
boundary conditions. The systems are initialized in a crystalline state and run for lOOr, at 
the frequencies / = lr“^, / = 3 t“^ and / = 5 t“^. The averages are calculated over 50 
and 250 strain cycles in the steady state respectively. Figure shows the parameters we 
have investigated in oscillatory shear conditions and the states we have identified. For a 
summary of the equilibrium state points under consideration the reader is referred to Fig. [T] 


The respective results are detailed in the Results section III 


B. Trajectory analysis 

a. Structure factors The static structure factor S{qx,qz) is evaluated in the velocity- 
vorticity plane at = 0. This reciprocal plane corresponds to a neutron experiment where 
the incident beam is parallel to the gradient direction of the shear flow. We calculate the 
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FIG. 3. (Color online) Out-of-equilibrium states [I: twinned FCC (a), II: disordered (o), III: 
sliding layers (■)] observed in the frequency-strain (/r-ymax) plane of the parameter space in our 
BD simulations at state point A, cf. Fig. The dashed lines are tentative boundaries of the 
respective nonequilibrium states. 


structure factor directly from the COM coordinates of the dumbbells as 

S{qo;,qz,0) = ( / ’ 


(14) 


, i=l 


where the angle brackets denote the trajectory average, which is taken in the steady state. 

b. Scattering intensity The scattering intensity of the suspension /(q) is the convolu¬ 
tion of the COM structure and the distribution of scattering centers within each particle. 
The scattering amplitude ^(q; u) describes the scattering of a single particle with orienta¬ 
tion u and constant internal density at the scattering vector q. Averaging the scattering 
amplitudes over all possible orientations yields the form factor P{q). The total scattering 
intensity is calculated as 


/(q) = (), 


( 15 ) 




where the scattering amplitude of a homogeneous dumbbell tilted by the angle 6 with respect 
to the scattering vector q is given by 


i 

A(q; u) = J dt cos (g cos 6R [t -|- L*]) 


-L* 


(16) 





c. Orientation In order to analyze the orientational behavior of the particles, a mea¬ 
sure of the orientation respecting the geometry of the system given is used. We plot the 
mean orientation with respect to the Cartesian axes, which coincide with the velocity, shear 
gradient and vorticity axes, as 

(^?)cyde (<) = {P2{cos0^))cck W, (V) 

where a denotes the coordinates x,y,z and Oa{t) the corresponding instantaneous angles. 
The average O^yde denotes sampling all time steps with the imposed strain state j{t). 

The orientational correlations with respect to the time are investigated in terms of di¬ 
rectional auto-correlation functions (DACFs) in the assumed steady state. Linearizing the 
equation of motion of the directors at some stationary state allows to define an effective 
rotational diffusion coefficient [l5]. The DACFs Ci{t) = (P/(u(t))), with I being the order 
of the Legendre polynomial Pi, may be approximated by 

Ci{t) = exp + 1);^^ • (18) 

d. Center of mass order parameters In order to distinguish crystalline structures, we 
use the local bond order analysis proposed by [10]. For each particle i a vector qfi) is 
defined by the components 

1 

= AT / - x ^ ^ (19) 

where y/TO(Rij) are the spherical harmonics for the normalized separation vectors Rjj, and 
Nh[i) is the number of the i-th particle’s neighbors. 
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where the second sum runs over all —I < mj < I fulfilling mi + m 2 + m 3 = 0 . 


9 




We define the ensemble averaged order parameters as 


(Qi) 

m 


and 


( 22 ) 

(23) 


where the angle brackets (...) denote the time average in the steady state. In this work 
the COM order is monitored in terms of averaged local order parameters {Q 4 ) and (114). 
These order parameters are sensitive to the configurations of the neighborhoods of solid-like 
particles. In particular, FCC and hexagonal close packed (HCP) structures are separated 
by a change of the sign of ( 114 ). 

e. Radial distribution function and Enskog collision rate The radial distribution func¬ 
tion (RDF) of the beads is defined as 


9(r) = 



(24) 


In a homogeneous and isotropic system the RDF depends only on the distance r = |r|. 

From a structural point of view, it is interesting to investigate the collision probability 
(or rate) during the shear and the phase transitions. It gives a rough picture of the average 
configurational freedom of the dumbbells under shear. The Enskog collision rate for hard 
spheres with diameter a is given by 071 


Fe = 9(cT)ro, 


(26) 


where g{a) is the contact value of the RDF, and Fq is the collision frequency in the dilute gas 
which is the ratio of the mean velocity and the free path in a suspension of hard spheres. In 
the Enskog approximation, which neglects correlated collisions, the self-diffusion coefficient 
is inverse proportional to the contact value 


De = 


2?7ir E 


(26) 


for a particle of mass m. Since our particles have smooth potentials, the RDE does not have 
its first maximum exactly at r = cr. Therefore, we take the height of the first maximum 
^(^"max) as a measure for the collision rate. 
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III. RESULTS AND DISCUSSION 


A. Steady-state structures and transitions 


We start with the discussion of BD simulated scattering intensities for selected state 
points. In Fig.|^the scattering intensities /(q) including the scattering amplitude calculated 
as defined in equations Eq. (15) and Eq. (jlhl) are shown in the velocity-vorticity plane. The 


data corresponds to state point B, cf. Fig. [1} at a frequency of 5r ^ for selected strains. 


The intensity plots illustrate the transition from the low strain twinned crystal state (I) 
to the high strain sliding layer state (III), corroborating with earlier results on a smaller 
aspect ratio In all ordered states a hexagonally ordered plane is parallel to the velocity- 
vorticity plane in the shear flow framework. A common feature of the scattering plots is the 
fact that the maxima may be seen on rings of constant scattering vector magnitude. At low 
strain amplitudes (Fig. |^,b) a shear twinned FCC crystal is stable, which has a densely 
packed direction perpendicular to the velocity direction, while at high strains (Fig. iJ-f)a 
dense direction parallel to the velocity is favorable. Here, the fully ordered systems tend 
to form two-dimensional HCP layers in the velocity-vorticity plane. These planes are the 
most densely packed crystallographic planes for FCC and HCP crystal structures and their 
stacking sequence determines the crystallographic type. 


Similar transitions had been observed in spherical systems [ZIIHIIII]. However, slightly 
anisotropic dumbbells introduce an additional, orientational degree of freedom, and, al¬ 
though weak, it leads to a more abrupt transition compared to the spherical case, as we 
have previously demonstrated for a smaller aspect ratio [37]. In equilibrium terms it could 
be said the mild anisotropy qualitatively change the transition from being continuous to 
discontinuous. This phenomenon is already present in the series of scattering intensities, 
as we observe a fully molten state for intermediate states, cf. Fig. B:. In contrast, for the 
hard spherical reference case, we always detect order in form of crystalline hybrids [H7] . 
At a strain amplitude about ~ 0.2 the system does not show long-range correlation 
and is nearly isotropic (Fig. |^). Close to this transition an anisotropy may be observed at 
very small scattering vectors in the reciprocal velocity-vorticity plane. Above this isotropic 
state hexagonal layers perform a zig-zag trajectory, which is indicated by a 30 ° tilt of the 
scattering pattern. In the present case we conclude from comparison to the idealized pic- 


11 





q^a (velocity) 



-10 -5 0 5 10 


Qxcr (velocity) 



FIG. 4. (Color online) Scattering intensities /(q) (250 cycles averaged) in the = 0 plane of 
sheared dumbbell suspensions (L* = 0.30,0 = 0.55, state point B) at frequency / = 5r“^. From 
top left to bottom right: increasing strain amplitude ymax^ 0.05 (a), 0.1 (b), 0.2 (c), 0.3 (d), 0.5 
(e) 1.0 (f) 


ture and from the reasoning concerning the volume fractions given by Ackerson [S] that 
the COM motion of the dumbbells follow strongly registered trajectories while maintaining 
two-dimensional hexagonal in-plane order. 


In the following, we investigate the particles’ orientational and translational structure 
under the action of shear. Figure [^summarizes the structural information in terms of aver¬ 
aged order parameters for state points A and B, cf. Fig. [T| at a frequency of 5 T“h Figurej^ 
shows that the average (Pf), quantifying the average dumbbell orientation along the flow 
direction, is slightly increasing with strain for both aspect ratios in the FCC state until it 
drops to zero in the disordered state and is non-zero again in the high-strain regime. In 
Fig.|§3 the translations order, described by {Qa), is shown for the dumbbells compared to 
the nearly hard sphere reference case (S). Thus, the non-monotonicity of the orientational 
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measure with respect to the strain amplitude is due to the loss of long-range order. The 
latter is not preserved in the transition region for sufficiently high aspect ratios, at which 
neighboring layers of different orientation cannot pass smoothly any more. The behavior at 
the transition shows clearly that the orientational and translational order changes abruptly 
in a discontinuous fashion for dumbbells, while, in contrast, the transition in suspension 
of hard spheres has a continuous character. This behavior is retained for dumbbells with 
a slender elongation of L* = 0.10 (system C), which show a ordered state at all investi¬ 
gated strain amplitudes. The average (P|^) for state point C shows a slight increase at low 
strain amplitudes as well approaching a plateau value of about 0.02 in a rather monotonous 
fashion beyond y^ax ^ 0.3. This corresponds well to our previous finding that long-time 
stress correlations due become important for dumbbells above approximately L* = 0.15 [18] . 
Additionally, the negative (WT) in state I indicates FCC dominated structure, in state III 
this parameter is vanishes on average indicating loss of FCC order. 

Moreover, the steady state structures and transitions depend on the shear frequency, 
which is summarized in the nonequilibrium phase diagram in Fig. The diagram shows 
the phases depicted for various strain amplitudes y^ax versus frequency / for state point A. 
In our simulations, state III does not appear at / = 1/r whereas state I is stable up to 
7 max ~ 0.30. For higher frequencies / = 3 /t and / = 5/r we observe the ordered state III 
above an amplitude of approximately 0.3, moving slightly to lower strains on increasing 
frequencies, which is sketched by the dashed lines in Fig. 

Let us now characterize the respective nonequilibrium states in more structural detail. 
Twinned-FCC regime (I): At low strain amplitudes a shear-twinned FCC dominated struc¬ 
ture is observed in the steady state. Figure]^ shows a series of snapshots taken at distinctive 
points in the strain cycle at which the instantaneous strain is minimal (y(t) = —ymax), zero 
(y(t) = 0) or maximal respectively. In this state, particles may follow the oscillatory flow 
and transfer between the triangular voids offered by neighboring layers. These void spaces 
are accessible in the vicinity of the extrema of the strain cycle, in between the particles are 
forced to pass particles of adjacent layers closely, which is referred to as bridge-stacking. 
This behavior finds expression in the transient orientation which is subtly coupled to the 
strain cycle through particle interaction. At these low strain amplitudes, y^ax < 0.1, the 
orientation shows a subtle interplay between the velocity and the vorticity axes while the 
amplitude in both is very small. This is exemplified in Fig. [^, where we show the time- 
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FIG. 5. (Color online) Order parameters for state points A and B at / = 5 r“^. (a) Averaged 
orientational order parameter in flow direction (Pf), and (b) translational order characterized by 
{Qa) (fllled symbols) and (IF 4 ) (open squares). 

resolved orientation in one shear cycle. In approaching the transition to the high strain 
state, 0.1 < 7max < 0.2, a slight separation of the velocity axis is observed, cf. Fig. [^, 
while the amplitudes of the {P2)cyc\e (-^Dcycie order of the average. 

The shear twinned FCC has a signature in the structure factor, see Fig. |^,b: the inner 
peaks are forbidden for equilibrium FCC crystals, in the present case we clearly observe 
non-vanishing peaks on the flrst ring and we see that their magnitude grows on increasing 
strain. Additionally, let it be noted, that the peaks on the velocity axis (g^ = 0) are the 
flrst to rise at very low strain amplitudes. Although the shear twin is not fully developed 
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at strains less than 0.2, the scattering intensity reveals that the equilibrium crystal is dis¬ 
turbed sufficiently from y^ax ^ 0.05. Let it be noted that the present very dense systems are 
crystalline in equilibrium and it has been shown that the crystalline state becomes stable 
for lower amplitudes on increasing density [9]. 



■» X 




FIG. 6 . (Color online) Simulation snapshots at prominent points in the shear cycle for 7 max = 
0.3, / = for state point A in the twinned FCC state I at (a) maximum, (b) zero, and (c) 

minimum instantaneous strain in the shear cycle. The particle radii are scaled by 1/2 for a better 
view of the structure. 



j_^_I_I_^_ II _I_I_I_^^_ II _I_ I I I ^ 

0 0.1 0.2 0 0.1 0.2 0 0.1 0.2 

tc/r tdr ted 


FIG. 7. (Color online) Time-resolved orientation within one cycle (averaged over 250 cycles) at 
frequency / = 5r“^, elongation L* — 0.30 (state point B), volume fraction </> = 0.55, and for 
different strain amplitudes: (a) 7 max = 0.05, (b) 7 max = 0.10, and (c) 7 max = 0.15. 

Intermediate disordered state (II): While for suspensions of hard spheres a disordered 
state in between the low-strain twinned-FCC (I) and high-strain sliding layer (III) regimes 


15 























is not observed, in fact, we find stable hybrid structures on reduction of the anisotropy, for 
sufiiciently elongated dumbbells the low-strain structure always melts fully at intermediate 
strain amplitudes. We find neither long-ranged translational order nor any orientational cor¬ 
relations whatsoever, which is confirmed considering the scattering intensities in Fig. &• A 
representative snapshot and cycle-averaged orientations are shown in Fig. clearly demon¬ 
strating disorder. Evidently, this is a distinctive behavior introduced by the orientational 
degree of freedom of the particles with a sufficient elongation. 



FIG. 8. (Color online) System (B) in the fnlly disordered state II at 7max = 0.2, snapshot (a) and 
cycle averaged orientational order parameters {P 2 ) cycle (^c) (b). 
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Sliding layer regime (III): While the centers of masses perform a zig-zag motion, a ten¬ 
dency of the particles’ orientation towards the velocity axis is observed, directly after reach¬ 
ing the critical strain to assemble into velocity oriented layers. Figure shows the cycle 
averaged orientation in the high strain regime, where a dense direction of each layer is 
aligned with the velocity direction. Here, the velocity and gradient directions of the direc¬ 
tors are clearly modulated by the shear cycle, while the vorticity direction is essentially flat. 
In this case, the modulus of the cycle average is maximal at times, when the instanta¬ 
neous strain vanishes (jit) = 0), and its amplitude is about 0.01 and the average 0.06. On 
average the directors slightly tend to be parallel to the velocity and perpendicular to the 
gradient. This coincides with a decoupling of the orientation from the imposed strain j{t), 
where the amplitudes of the (T’^)cycie significantly smaller than their respective averages. 


B. Kinetic properties 


Let us now turn our focus to rotational relaxation properties and particle collision rates. 
Figure [lO^ shows the inverse relaxation constant obtained from the exponential decays of 


the orientational autocorrelation functions according to Eq. (18) normalized by their short¬ 


time values at infinite dilution in absence of any external field. The effective rotational 
decay is enhanced in the disordered state while it is similar to the equilibrium case in the low 
strain regime and a significant slowing down is observed in the high strain state. At small 
strain amplitudes the values are slightly smaller than unity as we expect from the analysis of 
the equilibrium behavior with respect to volume fraction and elongation [IS] . In the regularly 
structured twinned FCC regime (I) the diffusion is basically constant with increasing strain 
amplitude. In the transition region 0.2 < y^ax < 0.3 we observe an elevated orientational 
diffusion where it is steeply curved with respect to the strain amplitude. On entering the 
high strain regime it jumps back to a value close to its initial value at rest. With increasing 
strain amplitude the orientational diffusion then slightly decreases. In the high strain regime 
we observe a state showing enhanced coupling of the dumbbells’ orientations in space and 
time. For the higher aspect ratio (system B) we observe similar behavior, while the peak in 
the disordered state (II) at about y^ax ~ 0.25 is much less pronounced. Also, on increasing 
the elongation, the normalized inverse time scales are smaller than in the former system (A). 
This can be explained considering the packing effect from our previous study [IH] , where we 
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FIG. 9. (Color online) Snapshot in the flow-gradient (x — y) plane (a) and oriental order parameters 
(b) in velocity (x), gradient (y) and vorticity (z) directions versus strain cycle at 7max = 1-00 for 
state point B in the high shear state III, where the orientations exhibit a flnite order modulated 
by the imposed shear. 

show that packing gets important for long-time relaxation from L* > 0.3. At the low aspect 
ratio of state point C, which does not show a fully disordered state, we observe essentially 
unhindered rotational relaxation and virtually no influence of the shear amplitude. 

The transition behavior of the rotational diffusion corresponds well to the contact value 
analysis, allowing us to connect structure and kinetics. Figure [T^ shows the values of the 
RDF g{r) at contact (r = rmax) of system A with L* = 0.24, = 0.55. On approaching the 
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melting strain at about 7max = 0.2 the contact value g(a) = r^/ro rises slightly. At the 
transition from the ordered shear twinned system (I) to the disordered state (II) the contact 
value shows a distinct jump of about 10 %. A smaller jump is observed at the transition 
to the ordered high strain regime where the contact value is less than in in the disordered 
state and reaches a plateau at about (/(r^ax) ~ 2.1. Following the inverse relationship 


of diffusion and contact value in the Enskog approximation Eq. (26), thus the diffusivity 
increases drastically in the transition region (II). 




Tmax 


FIG. 10. (Color online) (a) Orientational relaxation constants from the decays of the Ci correlation 
functions for state points A, B and C, and (b) contact values of the radial distribution function 
g{r) on increasing strain amplitude for state point A. 
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C. Comparison to experiments 
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FIG. 11. (Color online) Phase eqnilibrium (a) of hard dumbbells at L* « 0.3 from crystallisation 
experiments. The experimental phase diagram (b) (denoted by red squares (□)) is compared with 
the prediction of MC simulations (solid black line (•) [27]) for L* — 0.24 and L* — 0.30. 

In the following we introduce novel experiments on the rheology for the ’colloidal nitrogen’ 
case at L* = 0.3 (B). The details on the experimental and synthesis procedures may be 
obtained from our previous papers [HSl ETj. Comparison to the prediction of the phase 
diagram in Fig. 0) shows that the coexistence region (flnid/PC) shrinks on increasing 
elongation and the phase boundaries of the experimental systems fit very well. A glance 
at Fig. im also confirms that the fraction of crystalline sample in biphasic region is linear 
in concentration, and, thus, the equilibrium properties of the experimental correspondent 
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FIG. 12. (Color online) (a) Dependence of G' and G" on increasing 7max for hard dumbbells with 
L* Ri 0.30 in the plastic crystal phase under oscillatory shear at / = 1 Hz (circles) and / = 5Hz 
(squares). This measurement is performed with the default setting (100s/point), (b) Shear moduli 
versus strain amplitude 7max for the same system phase that is measured under oscillatory shear 
of / = IHz with 500s/point. The hlled symbols denote G', while open symbols represent G”. 


of system B are well-defined. The yielding behavior of the hard dumbbells with L* 


0.30 


in the plastic crystalline phase is displayed in Fig. [T^ . The experimental data indicate 
kinetic differences as the rheology of the system A is retained by either the frequency or 
the number of shear cycles per point: While the yielding behavior of hard dumbbells with 
L* ~ 0.30 in the plastic crystalline phase displayed in Fig. [T^ exhibits one event under 
oscillatory shear at / = 1 Hz there are two yielding events at / = 5 Hz. Compared with the 
former frequency, the number of applied shear cycles is increased by five times within the 
same measurement time for the experiment at / = 5 Hz. It is necessary to mention that the 
double yielding event is observed in the oscillatory shear field with / = 1 Hz as well when 
the measurement time is prolonged by five times as shown in Fig. [T^. Based on these three 
sets of experiments, it is concluded that the hard dumbbells with L* ~ 0.30 in the plastic 
crystalline phase can show the same double yielding behavior as the hard dumbbells with 
L* 0.24, but the former needs more or faster oscillations to achieve the steady state. The 


system of hard dumbbells with L* 0.30 (B) is closer to the glassy state predicted by Zhang 


and Schweizer than that with L* ~ 0.24 (A) at the same volume fraction of 0.6. Due to the 
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sample 

Dr/S-^ 

Pe^ 

(aligned) 

Pe^ 

(twinned 

fee) 

Pe^ 

(interme¬ 

diate) 

Pe^ 

(sliding 

layers) 

L* 

^ 0.24 

6.05 

0.24 

0.52 

0.62 

1.20 

L* 

^0.30 

17.70 

0.29 

0.75 

0.89 

1.62 


TABLE II. Comparison of critical Pe,. at the formation of vorticity alignment structure, the twinned 
FCC, the intermediate structure and the partially oriented sliding layers in the hard dumbbells 
with L* « 0.24 and L* « 0.30. is obtained from DLS and DDLS measurements. 

expected slowdown of the dynamics in the vicinity of the glass transition in the former case, 
stronger and longer oscillations are required to induce the same structural change as that 
with L* ~ 0.30. 

As detailed above, in the BD simulations at constant frequency (/ = 5r“^) (Fig. 
we observe the transition at nearly the same strain amplitude for systems A and B, while 
for both cases this transition clearly occurs at lower amplitudes. Moreover, rheo-SANS 
experiments have been carried out to investigate the underlying shear-induced structural 
evolution that corresponds to the double yielding behavior. For a direct compilation of 
scattering data and rheology see Fig. |A1| in the appendix. The rheo-SANS experiments 
clearly demonstrate that the longer dumbbells (B) undergo the same structural evolution 
to those with L* 0.24 (A), corresponding to the simulation results. With the increasing 
applied shear strains, fully crystallized longer dumbbells undergo the phase transition from 
twinned FCC (I) to the partially orientated sliding layer state (III), while being disordered in 
between (II). The main difference in the experimental realization is that the longer dumbbells 
need larger Peclet numbers Pe,. to induce the same nonequilibrium states as compared with 
hard dumbbells with L* 0.24, which we summarize in Table [III It should be noted that the 
simulated frequencies and, thus, the Peclet numbers are approximately five times as high as 
the experimental parameters. In our simulations, the transition from II to III, occurs nearly 
at the same values Pe^ = 4.56 and 4.64 for systems A and B, respectively. 


IV. CONCLUSIONS 

Following up on our previous study EH, we have here provided more evidence that 
the mild anisotropy of dumbbell-shaped particles leads to qualitative changes in the na- 
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ture of the nonequilibrium phase transitions of plastic crystals of spherical colloids under 
oscillatory shear, in particular that the continuous transition observed in spherical sys¬ 
tems transforms into a discontinuous one. The latter phenomenon must be attributed to 
rotational-translational couplings [ 3 T], absent in nearly spherical systems that apparently 
lead to dramatic changes in the structural and stress relaxation behavior. In fact, we recently 
showed by equilibrium BD simulations that plastic dumbbells crystals exhibit a dramatic 
increase of the linear shear response for high packing fractions above a critical aspect ratio 
of about 0.15 |1H]. With respect to the sequence of shear-induced states, the type and the 
dynamics of the equilibrium to twinned crystal transition remains an interesting issue, which 
may stimulate a future study. The observed strong transitions have substantial implications 
for rheology and the yielding behavior of anisotropic colloidal crystals, as detailed by |Chu 


the rheology, which we attribute to a dynamical slowing down. The present results also 
demonstrate that the thermosensitive dumbbell particles introduced before [HS] serve as an 
excellent and versatile model system for mildly anisotropic colloids to study their equilibrium 
and nonequilibrium structural and phase behavior. 


et al. 


The new experimental results also show a frequency and time dependency of 


Appendix A: Compilation of Rheology and Scattering 
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FIG. Al. Yielding behavior of the hard dumbbells with L* 0.30 (the DPM_b microgels) in the 
fully crystalline phase = 0.60) in the oscillatory shear of 5 Hz, and the corresponding 2D 
scattering patterns are measured by SANS. Along the dependence of G' and G" on various strains: 
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